Introduction to R
By: Thomas Marke, University of Innsbruck
2023-10
Kurzeinführung in die Datenanlyse mit R
Dies ist eine kurze Einführung in R für (Noch-)Nichtprogrammierer zum Einlesen, Analysieren und Darstellen exemplarischer Daten. Das Ziel ist nicht, Euch zu perfekten R ProgrammiererInnen zu machen, sondern durch Aufzeigen einiger praktischer und einfacher Funktionen die Verwendung von Programmiersprachen/Entwicklungsumgebungen schmackhaft zu machen und etwas “die Scheu” zu nehmen diese zu verwenden!
Neben der hier dargestellten Funktionalität könnt Ihr sehr einfach weitere Informationen zum Arbeiten mit R im Internet finden, z.B. durch Verwendung der Google Suchfunktion und den Ergebnissen auf der Seite stackoverflow.com - einfach mal ausprobieren.
Die hier verwendeten Beispieldaten sind Simulationen eines Energie-Bilanz-Schneemodels für den Standort Col de Porte (Frankreich) für die Saison 2005/2006 sowie Beobachtungen des Schneewasseräquivalents für den selben Standort und Zeitraum. Informationen zur Station und den Daten findet Ihr hier:
https://essd.copernicus.org/preprints/5/29/2012/essdd-5-29-2012.pdf
Die hier verwendeten Daten können hier heruntergeladen werden:
1 - Arbeitsumgebungen
Das Arbeiten mit R erfordert zunächst das R-Basispaket, welches für die jeweilig genutzten Betriebssysteme heruntergeladen und installiert werden muss. Pakete für die unterschiedlichen Betriebssysteme gibt es hier: https://cran.r-project.org/
Auf die dabei installierte Funktionalität (die R-Skriptsprache) kann man über unterschiedliche Umgebungen zugreifen:
R in der Konsole
Alle Betriebssysteme verfügen über eine “Konsole”, in dieser kann man verscheidene Kommandos eingeben um z.B. Verzeichnisse zu wechseln.
Befehl R: direktes Ausführen von R-Befehlen in der
Konsole
Befehl Rscript myscript.R: Ausführen bestehender
Skripte, hier “myscript.R”
R im Jupyter Notebook
Das Jupyter Notebook stellt eine übersichtliche Oberfläche zum Arbeiten in R dar, ist aber was die Datenverabreitung und den Funktionsumfang angeht in mancher Hinsicht etwas limitiert.
RStudio
RStudio ist die wohl am meisten verwendete Umgebung: https://www.rstudio.com/products/rstudio/download/
Hier kann man Code schreiben, Pakete installieren, Variablen bequem einsehen und Plots betrachten. Wir arbeiten in R-Studio, macht es doch gleich mal auf und öffnet ein neues Skript “Introduction.R” das Ihr im Ordner “R-Workspace” auf Eurem Rechner speichert.
Die Daten der Station “Col de Porte” (Frankreich) ladet Ihr auch gleich herunter und legt diese im gleichen Ordner ab!
2 - Grundlegendes
Programmablauf
Ein R-Script wird Zeile für Zeile von oben nach unten ausgeführt. So
können Variablen (z.B. Name) mit
Werten (z.B. Alex) gefüllt werden und für
numerische Variablen Berechnungen durchgeführt werden.
Kommentare
In jeder Programmiersprache ist es wichtig den Code
(die Befehlsabfolge) von Kommentaren (Erläuterungen die
ein Lesen des Codes vereinfachen und den Inhalt für Andere
nachvollziehbar machen) zu trennen. Kommentare können im R-Code durch
# eingeleitet werden - diese Teile werden dann nicht
ausgeführt, sondern sind nur für den Entwickler von Wert. Uns ist das
Kommentieren sehr wichtig, da es den LV-Leitern ermöglicht Euren Code zu
lesen und zu verstehen - bitte kommentiert doch deshalb Euren Code zur
besseren Nachvollziehbarkeit - das hilft auch dabei vor langer Zeit
geschriebene Skripte wieder verstehen/verwenden zu können.
Unten sehr Ihr ein Beispiel für einen Kommentar:
# Hello ausgeben
print('Hello')## [1] "Hello"Ausführen von Code
Code kann in R zeilenweise, für Bereiche oder für das gesamte Skript ausgeführt werden - hierfür gibt es fest definierte Shortcuts:
Ausführen der aktutellen Zeile/eines markierten Bereiches:
CTRL & Enter
Ausführen des gesamten Skriptes:
Shift & CTRL & s
Variablen
In der Programmierung arbeitet man mit “Variablen” die durch “Operationen” miteinander verarbeitet werden. So kann man einer frei zu definierenden Variablen einen Wert zuweisen und diese Variable dann mit einer anderen verrechnen.
Die Wertzuweisung erfolg in R durch
<- oder =
# Erste Berechnungen
apples = 11
oranges = 4
fruits = apples + orangesSchreiben wir den Namen der Variablen in eine Zeile und führen diese Zeile aus wird der Wert der Variable ausgegeben:
# Wert ausgeben
fruits## [1] 15Wir können uns den Typ einer Variablen durch
class() anzeigen lassen.
In diesem Fall ist es eine Variable vom Typ numeric,
also eine Zahl. Neben Zahlen gibt es zum Beispiel noch
character Variablen, diese enthalten eine Zeichenfolge.
“15” könnte also sowohl eine numerische Variable sein, aber auch eine
Zeichenfolge - mit Zeichenfolgen kann man aber natürlich nicht rechnen,
es empfiehlt sich bei Problemen im workflow immer mal den Typ einer
Variablen (v.a. bei Zahlen) zu checken, das erspart oft viel
Kopfzerbrechen:
# Klasse ausgeben
class(fruits)## [1] "numeric"Der Befehl str() gibt uns Information über die
Struktur unserer Daten.
In diesem Fall ist das nur die Information, das es sich um eine Zahl
mit dem Wert 15 handelt. Wir werden später mit DataFrames
arbeiten, hier können unterschiedliche Spalten unterschiedliche Typen
(z.B. Zahl, Zeichenfolge) haben.
# Struktur ausgeben
str(fruits)## num 15Einbinden von Bibliotheken (Libraries)
Während viele Befehle zur Standardfunktion von R gehören, erfordern einige Operationen Zusatzfunktionalität, welche in Form von Bibliotheken eingebunden werden kann.
Entsprechende Bibliotheken können in RStudio im Fenster rechts unten
durch Packages-->Install installiert werden, oder direkt
als R-Code durch den Befehl
install.packages("Paketname").
Wir installieren als Beispiel die Bibliothek this.path,
die uns später den Umgang mit Arbeitsverzeichnissen erleichtert
Um das Paket im Code verwenden zu können, muss der Befehl
library() in Zusammenhang mit dem Namen der einzubindenden
Bibliothek eingeben werden. Wir binden Testweise die Library
this.path ein:
library("this.path")In R arbeitet man oft mit vordefinierten Funktionen
die in Bibliotheken definiert sind. Will ich Hilfe zu einer solchen
Funktion, hilft Aufruf von help(Funktionsname). Wir testen
das anhand der Mittelwertsfunktion mean():
help(mean)Das Arbeitsverzeichnis
Habe ich irgendwo ein R-Skript angelegt, so ist es oft praktisch das Arbreitsverzeichnis auf den Ordner mit diesem Skript zu setzen. Da Daten in R standardmäßig (soweit nicht mit absolutem Pfad definiert) im Arbeitsverzeichnis gesucht werden, kann ich auf Daten die direkt im Arbeitsverzeichnis leigen ohne Angabe eines Pfads nur über den Dateinamen zugreifen. Liegen Daten in einem Unterverzeichnis des Arbeitsverzeichnisses, dann muss ich diesen Ordner relativ zum Arbeitsverzeichnis schon noch angeben.
Das Arbeitsverzeichnis kann durch den Befehl
setwd("/MeinOrdner") im R-Code gesetzt werden. Alternativ
kann in RStudio der aktuell angezeigte Ordner im Fenster rechts unten
als Arbeitsverzeichnis gesetzt werden
Files-->More-->Set As Working Directory.
Durch den Befehl getwd() kann ich mir das aktuelle
Arbeitsverzeichnis anzeigen lassen.
getwd()## [1] "/Users/t/pCloud Drive/Current/WS_2023/EU4 Gebirgsforschung/Introduction-R"Die eben installierte Bibliothek this.path ist in diesem
Zusammenhang besonders praktisch - sie ermöglicht es den Pfad in einem
R-Skript automatisch direkt auf den Pfad des Skriptes zu setzen. So muss
ich mir über die Lage meines R-Skriptes/Ordners auf meinem Rechner keine
Gedanken machen, solange alle benötigten Daten im gleichen Ordner
liegen.
# Bibliothek laden
library("this.path")
# Pfad des Skriptes auslesen und als Arbeitsverzeichnis setzen
script_path = this.path()
setwd(dirname(script_path))3 - Einlesen von Daten
Das Einlesen tabellarischer Daten in R ist sehr einfach. Durch die
Funktion read.table() können tabellarische Daten in
unterschiedlichen Formaten eingelesen werden. Man gibt hierfür den
Namen in ““, das Trennungszeichen
(z.B. sep=",") sowie die Information darüber, ob ein
Header (eine vorangehende Nichtdatenzeile mit
Spaltennamen) existiert an (z.B. header=TRUE). Weiter kann
ich angeben, welches Dezimaltrennzeichen in meinen
Daten verwendet wird (z.B. dec = "."). Mit dem Argument
skip= könnte man noch angeben, dass eine bestimmte Anzahl
an Zeilen übersprungen werden soll, z.B. wenn in den
ersten Zeilen unnötige Zusatzinformation enthalten ist.
Wir testen das mal mit einer Datei die Meteorologische Modellinputs
sowie eine Vielzahl von Modellautputs enthält. Die Datei sollte bereits
im Arbeitsverzeichnis liegen und heisst
ModelOutput.csv.
Tipp:
Schaut Euch tabellarischen Daten immer zunächst im Texteditor Eurer Wahl (Windows: Notepad, Mac: TextEdit, oder einfach in RStudio) an. So seht Ihr wie die Daten aussehen (Spalten, Header, Trennzeichen, etc.)!
Optionen für Trennzeichen sep="":
- Komma –>
"," - Strichpunkt –>
";" - Tabulator –>
"\t" - Kein Zeichen (white space) –>
""
# Daten einlesen
input_file = "ModelOutput.csv"
input_data = read.table(input_file, sep=",", dec = ".", header=TRUE)Datenstruktur
Die eingelesenen Daten befinden sich nun im Speicher des Rechners als
Werte der Variaben input_data und wir können uns diese
Daten einmal ansehen. Der Befehl names(input_data) zeigt
uns alle Spaltennamen an!
Tipp:
Unter Verwendung der Funktion
setNames(NameDataFrame,c("Spaltenname1","Spaltename2",...))
kann man die Spalten jederzeit umbenennen (machen wir später, hier aber
gerade nicht nötig).
# Spaltennamen ausgeben
names(input_data)## [1] "X" "Year" "Month" "Day" "Hour"
## [6] "ShortIn" "LongIn" "SolPrecip" "LiqPrecip" "AirTemp"
## [11] "RelHum" "WindSpeed" "SurfPress" "SWE" "SnowAlbedo"
## [16] "ShortOut" "LongOut" "SensFlux" "LatFlux" "ColdContent"
## [21] "Melt" "Runoff" "SnowDepth" "LiqWatCont" "IceCont"
## [26] "SnowTemp" "SnowDensity"Frage ich die Klasse der Variablen input_data ab, so
wird diese als data.frame angezeigt.
DataFrames sind sehr ähnlich einer Excel Tabelle mit
unterschiedlichen Zeilen und Spalten. Neben den DataFrames gibt es viele
andere Formate (z.B. das xts Format) die alle Ihre Vorteile und evtl.
leicht andere Befehlssyntax haben - aber bleiben wir bei unserem
DataFrame…
# Klasse ausgeben
class(input_data)## [1] "data.frame"Den gesamten Inhalt der Daten kann ich mir leicht ansehen, ich muss dazu nur alle Werte der Variablen ausgeben. Das mache ich indem ich die Variable ohne weiteren Befehl ausführe.
Wir könnten diese Zeile nach dem Anschauen unserer Daten im Code auch
mit # auskommentieren, damit der Output unseres Skriptes
nicht unnötig unübersichtlich wird :-)
# Alle Daten ausgeben
#input_dataWill ich mir nur die ersten 10 Zeilen
n=10 ansehen, kann ich die Funktion head()
einsetzen! Wir sehen hier einen Fehlwert (“Not Available” =
NA) in der Zeitreihe.
# Header (erste 10 Zeilen) ausgeben
head(input_data, n=10)## X Year Month Day Hour ShortIn LongIn SolPrecip LiqPrecip AirTemp RelHum
## 1 1 2005 10 1 0 0.0 283.1 0 0 277.8 78.2
## 2 2 2005 10 1 1 0.0 284.7 0 0 278.0 73.1
## 3 3 2005 10 1 2 0.0 285.8 0 0 277.7 76.1
## 4 4 2005 10 1 3 0.0 288.1 0 0 278.3 72.0
## 5 5 2005 10 1 4 0.0 293.9 0 0 277.7 73.0
## 6 6 2005 10 1 5 0.0 335.0 0 0 279.4 69.0
## 7 7 2005 10 1 6 1.4 347.2 0 0 280.7 57.8
## 8 8 2005 10 1 7 23.6 354.7 0 0 281.6 56.8
## 9 9 2005 10 1 8 96.4 346.0 0 0 282.9 59.8
## 10 10 2005 10 1 9 131.3 358.3 0 0 285.2 44.6
## WindSpeed SurfPress SWE SnowAlbedo ShortOut LongOut SensFlux LatFlux
## 1 0.6 87480 NA NA NA NA NA NA
## 2 0.0 87430 NA NA NA NA NA NA
## 3 1.0 87390 NA NA NA NA NA NA
## 4 0.5 87380 NA NA NA NA NA NA
## 5 0.2 87360 NA NA NA NA NA NA
## 6 1.2 87370 NA NA NA NA NA NA
## 7 1.1 87330 NA NA NA NA NA NA
## 8 1.1 87310 NA NA NA NA NA NA
## 9 0.6 87320 NA NA NA NA NA NA
## 10 0.6 87310 NA NA NA NA NA NA
## ColdContent Melt Runoff SnowDepth LiqWatCont IceCont SnowTemp SnowDensity
## 1 NA NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA NA
## 4 NA NA NA NA NA NA NA NA
## 5 NA NA NA NA NA NA NA NA
## 6 NA NA NA NA NA NA NA NA
## 7 NA NA NA NA NA NA NA NA
## 8 NA NA NA NA NA NA NA NA
## 9 NA NA NA NA NA NA NA NA
## 10 NA NA NA NA NA NA NA NAIch kann die einzelnen Spalten der Variable
ansprechen, indem ich das $ Zeichen verwende und
den Spaltennamen angebe! Wir testen das mit einem einfachen Ausdrucken
der Werte für die Stunde, wieder nur für die ersten 10 Zeilen!
# Header einer Spalte ausgeben
head(input_data$Hour, n=10)## [1] 0 1 2 3 4 5 6 7 8 9Datenstatistik
Oft interessiert uns die Statistik der Daten, z.B. der Mittelwert
oder die Standardabweichung - diese Information kann man sich über den
Befehl summary ausgeben lassen! Wir sehen wieder die
Information über enthaltene Fehlwerte NA in der
Datenreihe.
# Zusammenfassung ausgeben
summary(input_data)## X Year Month Day Hour
## Min. : 1 Min. :2005 Min. : 1.000 Min. : 1.00 Min. : 0.00
## 1st Qu.:1639 1st Qu.:2005 1st Qu.: 3.000 1st Qu.: 8.00 1st Qu.: 5.75
## Median :3276 Median :2006 Median : 5.000 Median :16.00 Median :11.50
## Mean :3276 Mean :2006 Mean : 6.033 Mean :15.68 Mean :11.50
## 3rd Qu.:4914 3rd Qu.:2006 3rd Qu.:10.000 3rd Qu.:23.00 3rd Qu.:17.25
## Max. :6552 Max. :2006 Max. :12.000 Max. :31.00 Max. :23.00
##
## ShortIn LongIn SolPrecip LiqPrecip
## Min. : 0.0 Min. :177.8 Min. :0.0000 Min. : 0.00000
## 1st Qu.: 0.0 1st Qu.:268.1 1st Qu.:0.0000 1st Qu.: 0.00000
## Median : 0.0 Median :297.4 Median :0.0000 Median : 0.00000
## Mean :102.3 Mean :292.0 Mean :0.0772 Mean : 0.05946
## 3rd Qu.:112.7 3rd Qu.:320.3 3rd Qu.:0.0000 3rd Qu.: 0.00000
## Max. :995.1 Max. :403.2 Max. :9.1080 Max. :10.26000
##
## AirTemp RelHum WindSpeed SurfPress
## Min. :258.3 Min. : 19.00 Min. :0.0000 Min. :84490
## 1st Qu.:270.1 1st Qu.: 64.80 1st Qu.:0.1000 1st Qu.:86460
## Median :275.7 Median : 81.30 Median :0.6000 Median :87000
## Mean :276.3 Mean : 76.86 Mean :0.9216 Mean :86856
## 3rd Qu.:281.7 3rd Qu.: 91.10 3rd Qu.:1.4000 3rd Qu.:87370
## Max. :297.0 Max. :102.20 Max. :8.5000 Max. :88410
##
## SWE SnowAlbedo ShortOut LongOut
## Min. : 0.02 Min. :0.5000 Min. :-741.89 Min. :-315.6
## 1st Qu.:148.58 1st Qu.:0.7385 1st Qu.: -60.54 1st Qu.:-314.1
## Median :287.95 Median :0.8586 Median : 0.00 Median :-297.9
## Mean :272.46 Mean :0.8240 Mean : -56.02 Mean :-295.5
## 3rd Qu.:392.83 3rd Qu.:0.9281 3rd Qu.: 0.00 3rd Qu.:-281.2
## Max. :499.71 Max. :0.9500 Max. : 0.00 Max. :-237.8
## NA's :2615 NA's :2615 NA's :2615 NA's :2615
## SensFlux LatFlux ColdContent Melt
## Min. :-67.595 Min. :-57.9717 Min. :-40.7333 Min. :0.0000
## 1st Qu.: 3.812 1st Qu.: -4.1091 1st Qu.:-11.0415 1st Qu.:0.0000
## Median : 13.787 Median : 0.0081 Median : -4.6945 Median :0.0000
## Mean : 16.071 Mean : 0.3702 Mean : -7.1116 Mean :0.1668
## 3rd Qu.: 26.453 3rd Qu.: 5.2245 3rd Qu.: -0.3866 3rd Qu.:0.0000
## Max. :119.472 Max. : 61.5705 Max. : 0.0000 Max. :4.6933
## NA's :2615 NA's :2615 NA's :2615 NA's :2615
## Runoff SnowDepth LiqWatCont IceCont
## Min. :0.0000 Min. :0.0002 Min. : 0.000 Min. : 0.0
## 1st Qu.:0.0000 1st Qu.:0.6069 1st Qu.: 0.000 1st Qu.:148.5
## Median :0.0000 Median :0.9808 Median : 0.000 Median :287.9
## Mean :0.1591 Mean :0.9455 Mean : 8.522 Mean :263.9
## 3rd Qu.:0.0000 3rd Qu.:1.2126 3rd Qu.:13.344 3rd Qu.:365.0
## Max. :5.0213 Max. :1.9417 Max. :49.973 Max. :488.0
## NA's :2615 NA's :2615 NA's :2615 NA's :2615
## SnowTemp SnowDensity
## Min. :254.5 Min. : 10.74
## 1st Qu.:265.4 1st Qu.:241.45
## Median :269.2 Median :270.24
## Mean :268.6 Mean :283.67
## 3rd Qu.:272.8 3rd Qu.:322.92
## Max. :273.1 Max. :473.33
## NA's :2615 NA's :26154 - Linien und Balkendiagramme
Einfache Diagramme (= Plots) zu erstellen ist in R relativ einfach.
Es gibt viele Plotting Bibliotheken (z.B. ggplot2), wir
wollen mit der Standard Plottingfunktion
plot() arbeiten, da diese für Einsteiger leichter zu
verstehen ist.
Für Fortgeschrittene ist das Plotten mit ggplot2 eine
tolle Alternative, da die Plots hier oft ohne viel Aufwand schöner
werden, im Internet gibt es hierzu tolle Tutorials, einfach mal googeln
:-)
Wir erstellen nun unseren ersten Plot und überlassen so gut wie alles
der plot() Funktion, d.h. wir geben ausser den X-Werten und
Y-Werten keine Argumente mit an die Funktion. Testweise plotten wir die
Lufttemperatur AirTemp gegen den Index X:
# Lufttemperatur plotten
plot(input_data$X, input_data$AirTemp)
Das funktioniert schonmal, allerdings hat dieser Plot noch einige Mängel die wir nun anpassen wollen.
Um unsere Zeitreihe auf der X-Achse gut darstellen zu können,
generieren wir nun eine Zeitinformation in unserem
DataFrame (wir nennen die neue Spalte
DateTime) aus den Spalten für Jahr, Monat, Tag und
Stunde. Um diese Zeitinformation zusammenzustellen verwenden
wir die R-Funktion paste() - sie ermöglicht das kombinieren
von Text und Variablenwerten.
Wichtig:
Variablen werden hier durch Komma getrennt angegeben, dabei ist
reiner Text in Anführungszeichen zu setzen, ganz am Ende wird nach einem
abschliessenden Komma mit sep='' definiert wie die
Einzeilteile zu trennen sind (hier ohne Inhalt weil wir keine
automatsiche Trennung, z.B. durch Leerzeichen wollen, diese fügen wir
lieber selbst ein).
# Spalte DateTime erzeugen
input_data$DateTime = paste(input_data$Year,'-',input_data$Month,'-',input_data$Day,' ',input_data$Hour,sep='')Wir sehen uns nun mit head() und str()
einmal an was in der Spalte DateTime abgelegt wurde:
# Header ausgeben
head(input_data$DateTime, n=10)## [1] "2005-10-1 0" "2005-10-1 1" "2005-10-1 2" "2005-10-1 3" "2005-10-1 4"
## [6] "2005-10-1 5" "2005-10-1 6" "2005-10-1 7" "2005-10-1 8" "2005-10-1 9"# Struktur ausgeben
str(input_data$DateTime)## chr [1:6552] "2005-10-1 0" "2005-10-1 1" "2005-10-1 2" "2005-10-1 3" ...Noch ist unsere Zeitinformation eine Zeichenfolge,
wir transformieren diese nun in eine POSIX
Zeitinformation, dabei geben wir das Format mit, in dem unsere
Zeichenfolge vorliegt, mit der Funktion head() geben wir
uns gleich die ersten Zeilen aus:
# In Zeitinformation transformieren und Header ausgeben
input_data$DateTime = as.POSIXct(input_data$DateTime,format="%Y-%m-%d %H")
head(input_data$DateTime, n=10)## [1] "2005-10-01 00:00:00 CEST" "2005-10-01 01:00:00 CEST"
## [3] "2005-10-01 02:00:00 CEST" "2005-10-01 03:00:00 CEST"
## [5] "2005-10-01 04:00:00 CEST" "2005-10-01 05:00:00 CEST"
## [7] "2005-10-01 06:00:00 CEST" "2005-10-01 07:00:00 CEST"
## [9] "2005-10-01 08:00:00 CEST" "2005-10-01 09:00:00 CEST"Nun plotten wir unsere Zeitreihe nocheinmal, diesmal die Lufttemperatur gegen unsere neue Zeitinformation:
# Lufttemperatur plotten
plot(input_data$DateTime, input_data$AirTemp)
Wir sehen, dass die plot-Funktion mit der Zeitinformation gut umgehen kann, ohne weitere Eingriffe werden in der X-Achse die Monate dargestellt.
Um weitere Mängel zu korrigieren plotten wir die Daten erneut -
diesmal geben der plot() Funktion noch ein paar
Informationen mehr mit, z.B.
dass wir nur Linien, keine Punkte darstellen wollen –>
type="l",dass wir gerne eine rote Linie hätten –>
col="red"unddass wir die Beschriftung der Achsen gerne selbst übernehmen –>
ann=FALSE.
# Lufttemperatur plotten
plot(input_data$DateTime, input_data$AirTemp, type="l", col="red", ann=FALSE)
Weitere Optionen für die Darsttellung der Daten
type="":
type=”l”–> nur Linientype=”p”–> nur Punktetype=”b”–> beidestype="h"–> vertikale Linien = Balkentype="n"–> nichts
Weiter Farboptionen col="":
redblueblack...
Einen schönen Überblick über die in R standardmäßig (es gibt viele Erweiterungspakete) verfügbaren Farben findet Ihr unter: http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
Nun wollen wir unseren Plot wie folgt erweitern:
als zweite Zeitreihe soll die simulierte Schneetemperatur
input_data$SnowTemphinzugefügt werdenein *Titel, Beschriftungen für die X- und Y-Achse sowie eine Legende soll hinzugefügt werden
wir wollen ein Grid hinter die Zeitreihen legen
Zweite Zeitreihe:
Wir fügen diese durch Aufruf von
lines()hinzu und verwenden das Attributltyum den Linetype (Linientyp) zu definieren. Wir geben der Lufttempertur den Stilsolid(-) und der Schneetemperatur den Stiltwodash(–). Beachtet, dass wir die Linientypen auch in der Legende angeben müssen, die Angaben erfolgen hier in Form eines Vektorsc()mit den beiden Linientypen. Als Farben wählen wirredundblue.Titel:
Wird eingefügt durch die Funktion
title(), sollte der Text zu lang sein kann man mit\neine neue Zeile einleiten!Legende:
Diese wird durch
legend()eingefügt, für die Legende gibt es verschiedene Platzierungsoptionen:"bottomright","bottom","bottomleft","left","topleft","top","topright","right"and"center"Durch Angabe eines Versatzes
inset=""in Prozent des Achsenbereichs sorgt man dafür, dass die Legende schön positioniert ist!Grid:
Hier verwenden wir die Funktion
grid(), der wir einen Linientyplty, eine Farbecolund eine Breitelwdmitgeben. Wir setzennxundnyaufNULL, wenn wir eine der beiden Achsen ohne Grid-Linien anzeigen wollten, würden wir diese aufNAsetzen.
Tipp:
Unterschiedliche Linientypen lty="":
twodashlongdashdotdashdotteddashedsolidblank
Diese können auch über Nummern angesprochen werden, wollt Ihr mehr über linetypes wissen besucht: http://www.sthda.com/english/wiki/line-types-in-r-lty
# Unser ursprünglicher Plot von oben, erweitert um die Angabe des Linientyps
plot(input_data$DateTime, input_data$AirTemp, type="l", col="red", ann=FALSE, lty='solid')
# Zweite Zeitreihe SnowTemp hinzufügen
lines(input_data$DateTime, input_data$SnowTemp, type="l", col="blue", lty='twodash')
# Diagrammtitel hinzufügen
title(main="Lufttemperatur und simulierte Schneetemperatur \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Monat")
title(ylab="Temperatur (K)")
# Grid hinzufügen
grid(lty = 2, col = "gray", lwd = 1)
# Legende hinzufügen, leichten Versatz angeben
legend("topleft", inset=.05, c("Lufttemperatur","Schneetemperatur"), col=c("red","blue"), lty=c('solid','twodash'))
Wir wollen nun noch den Niederschlag als
Balkendiagramm in blau plotten. Da wir am
Gesamtniederschlag interessiert sind und in unseren Daten fester
SolPrecip und flüssiger Niederschlag LiqPrecip
getrennt vorliegen, addieren wir diese zunächst in einer neuen Spalte
Precip. Dann plotten wir den Niederschlag unter Verwednung
des Darstellungstypes h (type="h").
# Flüssigen und festen Niederschlag in neuer Spalte zu Gesamtniederschlag addieren
input_data$Precip = input_data$LiqPrecip + input_data$SolPrecip
# Niederschlag plotten
plot(input_data$DateTime, input_data$Precip, type="h", col="blue", ann=FALSE)
# Diagrammtitel hinzufügen
title(main="Niederschlag \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Monat")
title(ylab="Stündlicher Niederschlag (mm)")
Als nächstes wollen wir einen Dichteplot erzeugen,
der uns zeigt, wie der Niederschlag über den Wertebereich verteilt ist.
Hierzu verwenden wir erneut die plot() Funktion, berechnen
uns aber mit der Funktion density() zunächst die
Dichteverteilung density_Precip, die wir einfach einmal
ausgeben um uns die Berechnungen einmal anzusehen:
# Dichtefunktion berechnen und ausgeben
density_Precip = density(input_data$Precip)
density_Precip##
## Call:
## density.default(x = input_data$Precip)
##
## Data: input_data$Precip (6552 obs.); Bandwidth 'bw' = 0.09585
##
## x y
## Min. :-0.2875 Min. :0.000000
## 1st Qu.: 2.4212 1st Qu.:0.000341
## Median : 5.1300 Median :0.001970
## Mean : 5.1300 Mean :0.092235
## 3rd Qu.: 7.8388 3rd Qu.:0.009271
## Max. :10.5475 Max. :3.677512Nun wollen wir uns die berechnete Dichtefunktion
density_Precip im Plot durch eine blaue Linie
(col="blue") darstellen, wieder wollen wir Titel und
Achsenbeschriftung selbst definieren (ann=FALSE).
# Plot erstellen
plot(density_Precip, col="blue", ann=FALSE)
# Diagrammtitel hinzufügen
title(main="Dichtefunktion des Niederschlags \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Niederschlag (mm)")
title(ylab="Dichte")
Wir fügen nun auch den festen Niederschlag hinzu, diesmal in einer
exotischeren Farbe, darkturquoise
# Ursprünglicher Plot für den festen Niederschlag
plot(density_Precip, col="blue", ann=FALSE)
# Festen Niederschlag hinzufügen, vorher Dichteverteilung berechnen
density_SolPrecip <- density(input_data$SolPrecip)
lines(density_SolPrecip, col="darkturquoise")
# Diagrammtitel hinzufügen
title(main="Dichtefunktion des Niederschlags \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Niederschlag (mm)")
title(ylab="Dichte")
Da ist die zweite Linie wohl etwas über das Ziel hinaus geschossen.
Das liegt daran, dass unser erster plot() Aufruf die
Y-Werte auf den kleineren Wertebereich optimiert hat. Also versuchen wir
es nochmal, diesmal begrenzen wir den Wertebereich der X- und
Y-Achse durch die Argumete xlim und
ylim. Die Angaben erfolgen hier in Form eines Vektors
c() mit den Minimal- und Maximalwerten. Grid und Legende
fügen wir auch gleich hinzu (beachtet, dass unser Grid zuerst
hinzugefügt wird, sonst liegt es über der Legende):
# Ursprünglicher Plot für den festen Niederschlag
plot(density_Precip, xlim=c(0,12), ylim=c(0,6), col="blue", ann=FALSE)
# Festen Niederschlag hinzufügen, vorher Dichteverteilung berechnen
density_SolPrecip <- density(input_data$SolPrecip)
lines(density_SolPrecip, col="darkturquoise")
# Diagrammtitel hinzufügen
title(main="Dichtefunktion des Niederschlags \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Niederschlag (mm)")
title(ylab="Dichte")
# Grid hinzufügen
grid(lty = 'dotted', col = "gray", lwd = 1)
# Legende hinzufügen
legend("topright", inset=.05, c("Niederschlag","Flüssiger Niederschlag"), col=c("blue","darkturquoise"), lty=c("solid","solid"))
5 - Subsets aus den Daten auswählen
Manchmal will man nur mit einem **Zeitbereich in *einer Datenreihe** arbeiten. Dann kann man auch Zeiträume extrahieren und am besten einer anderen Variablen zuordnen - so verlieren wir nicht unsere ursprünglichen Daten.
Das Auswählen von Daten aus einer Datenreiche kann über unterschiedliche Wege erfolgen:
Über den Index (z.B. Zeile 5-20)
Über eine Abfrage der Inhalte (z.B. Datum X - Y)
Über den Namen einer Spalte
Wir testen die Ansätze an einigen Beispielen unter Verwendung unseres
eingelesenen Datensatzes input_data.
Zunächst wollen wir alle Spalten des DataFrames, aber nur die
Zeilen 25 - 48 auswählen und in einen neuen DataFrame
index_subset schreiben, wir schauen diesen gleich mit
head() an:
# Bereich über den Index ausschneiden und Header ausgeben
index_subset = input_data[25:48,]
head(index_subset, n=50)## X Year Month Day Hour ShortIn LongIn SolPrecip LiqPrecip AirTemp RelHum
## 25 25 2005 10 2 0 0.0 339.7 0.000 0.4032 277.5 97.1
## 26 26 2005 10 2 1 0.0 322.8 0.000 0.0000 277.1 97.0
## 27 27 2005 10 2 2 0.0 331.4 0.000 0.0000 277.1 97.0
## 28 28 2005 10 2 3 0.0 327.2 0.000 0.0000 276.9 94.9
## 29 29 2005 10 2 4 0.0 326.7 0.000 0.0000 276.7 93.9
## 30 30 2005 10 2 5 0.0 328.6 0.000 0.7020 276.1 88.8
## 31 31 2005 10 2 6 0.3 325.0 0.000 0.9936 276.1 91.8
## 32 32 2005 10 2 7 3.9 328.2 0.000 3.1392 274.7 94.8
## 33 33 2005 10 2 8 14.4 317.7 0.000 2.4084 273.8 96.8
## 34 34 2005 10 2 9 35.3 324.4 0.000 1.4688 274.9 97.9
## 35 35 2005 10 2 10 53.1 327.2 0.000 3.4128 274.2 96.9
## 36 36 2005 10 2 11 43.3 314.4 4.248 0.0000 273.4 95.8
## 37 37 2005 10 2 12 122.8 318.5 0.000 2.6712 274.3 96.8
## 38 38 2005 10 2 13 104.4 325.9 0.000 2.5956 274.8 95.8
## 39 39 2005 10 2 14 63.3 322.8 0.000 3.1932 274.4 95.8
## 40 40 2005 10 2 15 46.4 324.3 0.000 2.4084 273.8 95.7
## 41 41 2005 10 2 16 25.8 317.9 0.000 3.2004 274.1 95.7
## 42 42 2005 10 2 17 12.2 323.0 0.000 2.8368 274.7 95.8
## 43 43 2005 10 2 18 0.3 326.4 0.000 1.1592 275.8 95.8
## 44 44 2005 10 2 19 0.0 328.6 0.000 1.5408 275.7 95.8
## 45 45 2005 10 2 20 0.0 329.7 0.000 2.5308 275.3 95.8
## 46 46 2005 10 2 21 0.0 328.6 0.000 0.6588 275.6 95.8
## 47 47 2005 10 2 22 0.0 330.0 0.000 0.0000 275.6 95.8
## 48 48 2005 10 2 23 0.0 330.6 0.000 0.2268 275.6 95.8
## WindSpeed SurfPress SWE SnowAlbedo ShortOut LongOut SensFlux
## 25 1.3 87070 NA NA NA NA NA
## 26 1.3 87050 NA NA NA NA NA
## 27 1.0 87010 NA NA NA NA NA
## 28 2.1 86980 NA NA NA NA NA
## 29 1.5 86950 NA NA NA NA NA
## 30 1.8 86960 NA NA NA NA NA
## 31 1.3 86950 NA NA NA NA NA
## 32 1.2 86970 NA NA NA NA NA
## 33 1.4 86980 NA NA NA NA NA
## 34 1.5 86970 NA NA NA NA NA
## 35 1.9 86990 NA NA NA NA NA
## 36 1.9 86970 4.2457159 0.9500000 -41.1350000 -315.637 1.725718
## 37 1.9 86920 4.3515649 0.9477556 -116.3843896 -315.637 7.938301
## 38 2.7 86910 3.8994986 0.9455224 -98.7125412 -315.637 13.828172
## 39 2.3 86890 3.5260823 0.9433004 -59.7109136 -315.637 9.552238
## 40 1.9 86890 3.2737728 0.9410894 -43.6665483 -315.637 4.486866
## 41 1.7 86890 3.0643143 0.9388895 -24.2233481 -315.637 6.206740
## 42 2.5 86890 2.6880580 0.9367005 -11.4277460 -315.637 12.417438
## 43 4.1 86890 1.9865034 0.9345224 -0.2803567 -315.637 29.062365
## 44 3.1 86910 1.3159252 0.9323552 0.0000000 -315.637 23.255057
## 45 3.5 86910 0.6610287 0.9301989 0.0000000 -315.637 21.195883
## 46 4.1 86890 NA NA NA NA NA
## 47 4.6 86900 NA NA NA NA NA
## 48 3.4 86910 NA NA NA NA NA
## LatFlux ColdContent Melt Runoff SnowDepth LiqWatCont IceCont
## 25 NA NA NA NA NA NA NA
## 26 NA NA NA NA NA NA NA
## 27 NA NA NA NA NA NA NA
## 28 NA NA NA NA NA NA NA
## 29 NA NA NA NA NA NA NA
## 30 NA NA NA NA NA NA NA
## 31 NA NA NA NA NA NA NA
## 32 NA NA NA NA NA NA NA
## 33 NA NA NA NA NA NA NA
## 34 NA NA NA NA NA NA NA
## 35 NA NA NA NA NA NA NA
## 36 -1.7987422 0 0.03744966 0.000000 0.041648549 0.03744966 4.2082662
## 37 3.8052949 0 0.28658825 2.570183 0.041871428 0.42505480 3.9265101
## 38 7.0398523 0 0.47200139 3.056606 0.036821677 0.43605044 3.4634482
## 39 3.9392064 0 0.33281813 3.571618 0.032688916 0.39045008 3.1356322
## 40 0.2261266 0 0.21478352 2.660997 0.029809376 0.35263695 2.9211359
## 41 1.7381866 0 0.18662678 3.412066 0.027416347 0.32759800 2.7367163
## 42 6.0733456 0 0.36357316 3.220769 0.023640441 0.30720265 2.3808553
## 43 18.5137561 0 0.68901821 1.884264 0.017179446 0.27115675 1.7153466
## 44 14.5996749 0 0.61846496 2.229917 0.011194638 0.20050427 1.1154210
## 45 12.3960564 0 0.60329986 3.201438 0.005533609 0.13316663 0.5278621
## 46 NA NA NA NA NA NA NA
## 47 NA NA NA NA NA NA NA
## 48 NA NA NA NA NA NA NA
## SnowTemp SnowDensity DateTime Precip
## 25 NA NA 2005-10-02 00:00:00 0.4032
## 26 NA NA 2005-10-02 01:00:00 0.0000
## 27 NA NA 2005-10-02 02:00:00 0.0000
## 28 NA NA 2005-10-02 03:00:00 0.0000
## 29 NA NA 2005-10-02 04:00:00 0.0000
## 30 NA NA 2005-10-02 05:00:00 0.7020
## 31 NA NA 2005-10-02 06:00:00 0.9936
## 32 NA NA 2005-10-02 07:00:00 3.1392
## 33 NA NA 2005-10-02 08:00:00 2.4084
## 34 NA NA 2005-10-02 09:00:00 1.4688
## 35 NA NA 2005-10-02 10:00:00 3.4128
## 36 273.15 101.9415 2005-10-02 11:00:00 4.2480
## 37 273.15 103.9268 2005-10-02 12:00:00 2.6712
## 38 273.15 105.9023 2005-10-02 13:00:00 2.5956
## 39 273.15 107.8678 2005-10-02 14:00:00 3.1932
## 40 273.15 109.8236 2005-10-02 15:00:00 2.4084
## 41 273.15 111.7696 2005-10-02 16:00:00 3.2004
## 42 273.15 113.7059 2005-10-02 17:00:00 2.8368
## 43 273.15 115.6326 2005-10-02 18:00:00 1.1592
## 44 273.15 117.5496 2005-10-02 19:00:00 1.5408
## 45 273.15 119.4571 2005-10-02 20:00:00 2.5308
## 46 NA NA 2005-10-02 21:00:00 0.6588
## 47 NA NA 2005-10-02 22:00:00 0.0000
## 48 NA NA 2005-10-02 23:00:00 0.2268Das hat gut geklappt, aber was bedeutet die fehlende Angabe nach dem Komma?
Unser DataFrame ist in Zeilen und Spalten organisiert, die erste
Dimension beinhaltet die Zeilen, hier haben wir mit 25:48
einen Bereich gewählt, die zweite Dimension die Spalten, hier haben wir
keine Angabe gemacht, somit wurden alle Spalten gewählt!
Im Beispiel haben wir alle Stunden des 02.10.2005 ausgeschnitten! Aber natürlich muss ich hierzu vorher den zugehörigen Indexbereich kennen.
Um unabhängig von den Indizes ein Zeitfenster auszuwählen kann ich
auch den Zeitraum definieren. Das Datum setzen uns
wieder mit der paste-Funktion aus year, month
und day zusammen und transformieren das Ergebnis dann
gleich wieder in ein Datum:
# Datumsspalte erzeugen
input_data$Date = paste(input_data$Year,'-',input_data$Month,'-',input_data$Day,sep='')
input_data$Date = as.Date(input_data$Date)
# Bereich über Datum ausschneiden und Header ausgeben
time_subset = input_data[input_data$Date>"2005-10-01" & input_data$Date<"2005-10-03",]
head(time_subset, n=50)## X Year Month Day Hour ShortIn LongIn SolPrecip LiqPrecip AirTemp RelHum
## 25 25 2005 10 2 0 0.0 339.7 0.000 0.4032 277.5 97.1
## 26 26 2005 10 2 1 0.0 322.8 0.000 0.0000 277.1 97.0
## 27 27 2005 10 2 2 0.0 331.4 0.000 0.0000 277.1 97.0
## 28 28 2005 10 2 3 0.0 327.2 0.000 0.0000 276.9 94.9
## 29 29 2005 10 2 4 0.0 326.7 0.000 0.0000 276.7 93.9
## 30 30 2005 10 2 5 0.0 328.6 0.000 0.7020 276.1 88.8
## 31 31 2005 10 2 6 0.3 325.0 0.000 0.9936 276.1 91.8
## 32 32 2005 10 2 7 3.9 328.2 0.000 3.1392 274.7 94.8
## 33 33 2005 10 2 8 14.4 317.7 0.000 2.4084 273.8 96.8
## 34 34 2005 10 2 9 35.3 324.4 0.000 1.4688 274.9 97.9
## 35 35 2005 10 2 10 53.1 327.2 0.000 3.4128 274.2 96.9
## 36 36 2005 10 2 11 43.3 314.4 4.248 0.0000 273.4 95.8
## 37 37 2005 10 2 12 122.8 318.5 0.000 2.6712 274.3 96.8
## 38 38 2005 10 2 13 104.4 325.9 0.000 2.5956 274.8 95.8
## 39 39 2005 10 2 14 63.3 322.8 0.000 3.1932 274.4 95.8
## 40 40 2005 10 2 15 46.4 324.3 0.000 2.4084 273.8 95.7
## 41 41 2005 10 2 16 25.8 317.9 0.000 3.2004 274.1 95.7
## 42 42 2005 10 2 17 12.2 323.0 0.000 2.8368 274.7 95.8
## 43 43 2005 10 2 18 0.3 326.4 0.000 1.1592 275.8 95.8
## 44 44 2005 10 2 19 0.0 328.6 0.000 1.5408 275.7 95.8
## 45 45 2005 10 2 20 0.0 329.7 0.000 2.5308 275.3 95.8
## 46 46 2005 10 2 21 0.0 328.6 0.000 0.6588 275.6 95.8
## 47 47 2005 10 2 22 0.0 330.0 0.000 0.0000 275.6 95.8
## 48 48 2005 10 2 23 0.0 330.6 0.000 0.2268 275.6 95.8
## WindSpeed SurfPress SWE SnowAlbedo ShortOut LongOut SensFlux
## 25 1.3 87070 NA NA NA NA NA
## 26 1.3 87050 NA NA NA NA NA
## 27 1.0 87010 NA NA NA NA NA
## 28 2.1 86980 NA NA NA NA NA
## 29 1.5 86950 NA NA NA NA NA
## 30 1.8 86960 NA NA NA NA NA
## 31 1.3 86950 NA NA NA NA NA
## 32 1.2 86970 NA NA NA NA NA
## 33 1.4 86980 NA NA NA NA NA
## 34 1.5 86970 NA NA NA NA NA
## 35 1.9 86990 NA NA NA NA NA
## 36 1.9 86970 4.2457159 0.9500000 -41.1350000 -315.637 1.725718
## 37 1.9 86920 4.3515649 0.9477556 -116.3843896 -315.637 7.938301
## 38 2.7 86910 3.8994986 0.9455224 -98.7125412 -315.637 13.828172
## 39 2.3 86890 3.5260823 0.9433004 -59.7109136 -315.637 9.552238
## 40 1.9 86890 3.2737728 0.9410894 -43.6665483 -315.637 4.486866
## 41 1.7 86890 3.0643143 0.9388895 -24.2233481 -315.637 6.206740
## 42 2.5 86890 2.6880580 0.9367005 -11.4277460 -315.637 12.417438
## 43 4.1 86890 1.9865034 0.9345224 -0.2803567 -315.637 29.062365
## 44 3.1 86910 1.3159252 0.9323552 0.0000000 -315.637 23.255057
## 45 3.5 86910 0.6610287 0.9301989 0.0000000 -315.637 21.195883
## 46 4.1 86890 NA NA NA NA NA
## 47 4.6 86900 NA NA NA NA NA
## 48 3.4 86910 NA NA NA NA NA
## LatFlux ColdContent Melt Runoff SnowDepth LiqWatCont IceCont
## 25 NA NA NA NA NA NA NA
## 26 NA NA NA NA NA NA NA
## 27 NA NA NA NA NA NA NA
## 28 NA NA NA NA NA NA NA
## 29 NA NA NA NA NA NA NA
## 30 NA NA NA NA NA NA NA
## 31 NA NA NA NA NA NA NA
## 32 NA NA NA NA NA NA NA
## 33 NA NA NA NA NA NA NA
## 34 NA NA NA NA NA NA NA
## 35 NA NA NA NA NA NA NA
## 36 -1.7987422 0 0.03744966 0.000000 0.041648549 0.03744966 4.2082662
## 37 3.8052949 0 0.28658825 2.570183 0.041871428 0.42505480 3.9265101
## 38 7.0398523 0 0.47200139 3.056606 0.036821677 0.43605044 3.4634482
## 39 3.9392064 0 0.33281813 3.571618 0.032688916 0.39045008 3.1356322
## 40 0.2261266 0 0.21478352 2.660997 0.029809376 0.35263695 2.9211359
## 41 1.7381866 0 0.18662678 3.412066 0.027416347 0.32759800 2.7367163
## 42 6.0733456 0 0.36357316 3.220769 0.023640441 0.30720265 2.3808553
## 43 18.5137561 0 0.68901821 1.884264 0.017179446 0.27115675 1.7153466
## 44 14.5996749 0 0.61846496 2.229917 0.011194638 0.20050427 1.1154210
## 45 12.3960564 0 0.60329986 3.201438 0.005533609 0.13316663 0.5278621
## 46 NA NA NA NA NA NA NA
## 47 NA NA NA NA NA NA NA
## 48 NA NA NA NA NA NA NA
## SnowTemp SnowDensity DateTime Precip Date
## 25 NA NA 2005-10-02 00:00:00 0.4032 2005-10-02
## 26 NA NA 2005-10-02 01:00:00 0.0000 2005-10-02
## 27 NA NA 2005-10-02 02:00:00 0.0000 2005-10-02
## 28 NA NA 2005-10-02 03:00:00 0.0000 2005-10-02
## 29 NA NA 2005-10-02 04:00:00 0.0000 2005-10-02
## 30 NA NA 2005-10-02 05:00:00 0.7020 2005-10-02
## 31 NA NA 2005-10-02 06:00:00 0.9936 2005-10-02
## 32 NA NA 2005-10-02 07:00:00 3.1392 2005-10-02
## 33 NA NA 2005-10-02 08:00:00 2.4084 2005-10-02
## 34 NA NA 2005-10-02 09:00:00 1.4688 2005-10-02
## 35 NA NA 2005-10-02 10:00:00 3.4128 2005-10-02
## 36 273.15 101.9415 2005-10-02 11:00:00 4.2480 2005-10-02
## 37 273.15 103.9268 2005-10-02 12:00:00 2.6712 2005-10-02
## 38 273.15 105.9023 2005-10-02 13:00:00 2.5956 2005-10-02
## 39 273.15 107.8678 2005-10-02 14:00:00 3.1932 2005-10-02
## 40 273.15 109.8236 2005-10-02 15:00:00 2.4084 2005-10-02
## 41 273.15 111.7696 2005-10-02 16:00:00 3.2004 2005-10-02
## 42 273.15 113.7059 2005-10-02 17:00:00 2.8368 2005-10-02
## 43 273.15 115.6326 2005-10-02 18:00:00 1.1592 2005-10-02
## 44 273.15 117.5496 2005-10-02 19:00:00 1.5408 2005-10-02
## 45 273.15 119.4571 2005-10-02 20:00:00 2.5308 2005-10-02
## 46 NA NA 2005-10-02 21:00:00 0.6588 2005-10-02
## 47 NA NA 2005-10-02 22:00:00 0.0000 2005-10-02
## 48 NA NA 2005-10-02 23:00:00 0.2268 2005-10-02Um nun bestimmte Variablen aus dem
DataFrame auszuwählen, kann ich diese bei der Auswahl wie
folgt ansprechen (hier am Beispiel des Schneewasseräquivalents
SWE), unsere Zeitinformation DateTime nehmen
wir gleich mit:
# Bereich über Spaltennamen auswählen
variable_subset = input_data[,c("DateTime","SWE")]
head(variable_subset, n=50)## DateTime SWE
## 1 2005-10-01 00:00:00 NA
## 2 2005-10-01 01:00:00 NA
## 3 2005-10-01 02:00:00 NA
## 4 2005-10-01 03:00:00 NA
## 5 2005-10-01 04:00:00 NA
## 6 2005-10-01 05:00:00 NA
## 7 2005-10-01 06:00:00 NA
## 8 2005-10-01 07:00:00 NA
## 9 2005-10-01 08:00:00 NA
## 10 2005-10-01 09:00:00 NA
## 11 2005-10-01 10:00:00 NA
## 12 2005-10-01 11:00:00 NA
## 13 2005-10-01 12:00:00 NA
## 14 2005-10-01 13:00:00 NA
## 15 2005-10-01 14:00:00 NA
## 16 2005-10-01 15:00:00 NA
## 17 2005-10-01 16:00:00 NA
## 18 2005-10-01 17:00:00 NA
## 19 2005-10-01 18:00:00 NA
## 20 2005-10-01 19:00:00 NA
## 21 2005-10-01 20:00:00 NA
## 22 2005-10-01 21:00:00 NA
## 23 2005-10-01 22:00:00 NA
## 24 2005-10-01 23:00:00 NA
## 25 2005-10-02 00:00:00 NA
## 26 2005-10-02 01:00:00 NA
## 27 2005-10-02 02:00:00 NA
## 28 2005-10-02 03:00:00 NA
## 29 2005-10-02 04:00:00 NA
## 30 2005-10-02 05:00:00 NA
## 31 2005-10-02 06:00:00 NA
## 32 2005-10-02 07:00:00 NA
## 33 2005-10-02 08:00:00 NA
## 34 2005-10-02 09:00:00 NA
## 35 2005-10-02 10:00:00 NA
## 36 2005-10-02 11:00:00 4.2457159
## 37 2005-10-02 12:00:00 4.3515649
## 38 2005-10-02 13:00:00 3.8994986
## 39 2005-10-02 14:00:00 3.5260823
## 40 2005-10-02 15:00:00 3.2737728
## 41 2005-10-02 16:00:00 3.0643143
## 42 2005-10-02 17:00:00 2.6880580
## 43 2005-10-02 18:00:00 1.9865034
## 44 2005-10-02 19:00:00 1.3159252
## 45 2005-10-02 20:00:00 0.6610287
## 46 2005-10-02 21:00:00 NA
## 47 2005-10-02 22:00:00 NA
## 48 2005-10-02 23:00:00 NA
## 49 2005-10-03 00:00:00 NA
## 50 2005-10-03 01:00:00 NA6 - Daten aggregieren
Oft liegen Daten zeitlich feiner aufgelöst (z.B. als Stundenwerte) vor als wir sie bräuchten (z.B. Monatswerte). Um die Daten zeitlich zu aggregieren, braucht es ein Zeitkriterium anhand dessen aggregiert werden kann, z.B. eine Spalte die nur Jahre- und Monatsinformation enthält - ich kann dann z.B. den Mittelwert oder die Summe über dieses Kriterium berechnen.
Wir testen das anhand unserer Niederschlagsdaten, dafür fügen wir
unserem DataFrame eine Spalte hinzu, die nur Jahr und Monat
enthält und verwenden die Funktion yearmon() des Paketes
zoo um diese Zeichenfolge in eine Zeitinformation zu
transformieren:
# Bibliothek einbinden
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# Spalte mit Jahr und Monat erstellen
input_data$YearMonth = paste(input_data$Year,'-',input_data$Month,sep='')
input_data$YearMonth = as.yearmon(input_data$YearMonth,format="%Y-%m")
# Header ausgeben
head(input_data$YearMonth, n=10)## [1] "Oct 2005" "Oct 2005" "Oct 2005" "Oct 2005" "Oct 2005" "Oct 2005"
## [7] "Oct 2005" "Oct 2005" "Oct 2005" "Oct 2005"Nun verwenden wir die aggregate() Funktion um die
Niederschlagssumme für die Monate zu berechnen, das
Ergebnis schauen wir gleich mit head() an:
# Monatliche Niederschläge berechnen
monthly_Precip <- aggregate(input_data$Precip,FUN = sum,by = list(Group.date = input_data$YearMonth))
# Header ausgeben
head(monthly_Precip, n=10)## Group.date x
## 1 Oct 2005 164.82636
## 2 Nov 2005 99.27468
## 3 Dec 2005 158.77692
## 4 Jan 2006 97.14528
## 5 Feb 2006 109.91376
## 6 Mar 2006 136.73599
## 7 Apr 2006 26.43555
## 8 May 2006 49.80496
## 9 Jun 2006 52.51841Die Spaltenüberschriften passen wir noch so an, dass
man auch weiss was hier enthalten ist! Wir verwenden hierzu die Funktion
setNames():
# Spaltennamen definieren
monthly_Precip = setNames(monthly_Precip,c("Month","Precip"))
# Header ausgeben
head(monthly_Precip, n=10)## Month Precip
## 1 Oct 2005 164.82636
## 2 Nov 2005 99.27468
## 3 Dec 2005 158.77692
## 4 Jan 2006 97.14528
## 5 Feb 2006 109.91376
## 6 Mar 2006 136.73599
## 7 Apr 2006 26.43555
## 8 May 2006 49.80496
## 9 Jun 2006 52.51841Die hier erhaltenen Monatswerte stellen wir nun als Balkendiagramm in
blau dar. Verwendet dabei den Typ type="h" sowie das
Argument lwd um die Dicke Eurer Balken anzupassen:
Tipp:
Die Farbe kann wie gewohnt einfach unter Angabe einer Farbe gesetz
werden (col="blue") oder über einen RGB Wert definiert
werden:
col=rgb(0,0,0,alpha=0.3)
Diese Option bietet auch die Möglichkeit einen Transparenzwert
alpha zu definieren, allerdings ist aufzupassen, da RGB
Werte in der Regel einen Bereich von 0-255 haben, hier aber eine
Intensität von 0-1 eingesetzt werden muss! Man muss die RGB-Werte falls
im 0-255er Wertebereich vorliegend, also normieren, d.h. durch den
Maximalwert teilen:
Für die Farbe blau heisst das: rgb(0,0,255) –>
rgb(0/255,0/255,255/255)
Eine Übersicht über verschiedene Farben und deren RGB Werte findet Ihr unter http://www.farb-tabelle.de/de/farbtabelle.htm!
# Monatlichen Niederschlag plotten
plot(monthly_Precip$Month, monthly_Precip$Precip, type="h", col=rgb(0/255,0/255,255/255, alpha=0.5), ann=FALSE, lwd=20)
# Diagrammtitel hinzufügen
title(main="Monatlicher Niederschlag \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Monat")
title(ylab="Niederschlag (mm)")
Der Plot sieht soweit ganz gut aus, lediglich die Zeitachse ist noch
eine Mischung aus Monat und Jahr - R tut das um uns den Jahreswechsel zu
zeigen, aber das sieht nicht so schön aus. Wir werden die Formatierung
der Zeitachse also selbst in die Hand nehmen, dafür geben wir in unserem
plot-Aufruf den Parameter xaxt="n" mit so dass zunächst
keine Achsenbeschriftung erfolgt, dann fügen wir am Schluss ein eigenes
Axen-Element mit entsprechend formatierter Zeitinformation über den
Befehl axis() hinzu:
# Monatlichen Niederschlag plotten
plot(monthly_Precip$Month, monthly_Precip$Precip, type="h", col=rgb(0/255,0/255,255/255, alpha=0.5), xaxt = "n", ann=FALSE, lwd=20)
# Diagrammtitel hinzufügen
title(main="Monatlicher Niederschlag \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Monat")
title(ylab="Niederschlag (mm)")
# Add dates to x-axis
axis(1,monthly_Precip$Month,format(monthly_Precip$Month, "%Y-%m"))
7 - Scatterplots
Scatterplots (deutsch: Streudiagramme) ermöglichen es den Zusammenhang zweier Größen darzustellen. Dabei wird die abhängige Variable (Y) als Funktion der unabhängigen Variable (X) dargestellt.
Wir wollen nun Scatterplots auf Basis simulierter und beobachteter Schneegrößen erstellen - diese Art der Auswertung stellt eine wichtige Säule in der Evaluation von Modellergebnissen dar.
Wir lesen dazu beobachte Werte des Schneewasseräquivalents zum
Vergleich mit den Simulationen ein - diese Daten befinden sich in der
Datei Observations.txt im Arbeitsverzeichnis.
Schaut Euch die Daten zunächst im Editor an und versucht den Lesebefehl richtig zu konfigurieren, ohne unten nachzusehen!
# Beobachtungsdaten einlesen
observation_file = "Observations.txt"
daily_observed_SWE = read.table(observation_file, sep="", dec = ".", header=TRUE)Wir sehen uns das Ergebnis des Einlesens mit head()
an:
# Header ausgeben
head(daily_observed_SWE, n=10)## year month day SWE
## 1 2005 10 1 0
## 2 2005 10 2 0
## 3 2005 10 3 0
## 4 2005 10 4 0
## 5 2005 10 5 0
## 6 2005 10 6 0
## 7 2005 10 7 0
## 8 2005 10 8 0
## 9 2005 10 9 0
## 10 2005 10 10 0Wir sehen dass alles geklappt hat, allerdings liegen die
Beobachtungen als Tageswerte vor - wir müssen also zunächts die
Simulationen auf Tageswerte aggregieren. Dafür
verwenden wir wieder die aggregate() Funktion zusammen mit
unserem vorhin erzeugten Datum:
# Stündliche SWE Simulationen auf Tageswerte aggregieren
daily_simulated_SWE <- aggregate(input_data$SWE,FUN = mean,by = list(Group.date = input_data$Date))Ein Blick auf das Ergebnis:
# Header ausgeben
head(daily_simulated_SWE, n=60)## Group.date x
## 1 2005-10-01 NA
## 2 2005-10-02 NA
## 3 2005-10-03 NA
## 4 2005-10-04 NA
## 5 2005-10-05 NA
## 6 2005-10-06 NA
## 7 2005-10-07 NA
## 8 2005-10-08 NA
## 9 2005-10-09 NA
## 10 2005-10-10 NA
## 11 2005-10-11 NA
## 12 2005-10-12 NA
## 13 2005-10-13 NA
## 14 2005-10-14 NA
## 15 2005-10-15 NA
## 16 2005-10-16 NA
## 17 2005-10-17 NA
## 18 2005-10-18 NA
## 19 2005-10-19 NA
## 20 2005-10-20 NA
## 21 2005-10-21 NA
## 22 2005-10-22 NA
## 23 2005-10-23 NA
## 24 2005-10-24 NA
## 25 2005-10-25 NA
## 26 2005-10-26 NA
## 27 2005-10-27 NA
## 28 2005-10-28 NA
## 29 2005-10-29 NA
## 30 2005-10-30 NA
## 31 2005-10-31 NA
## 32 2005-11-01 NA
## 33 2005-11-02 NA
## 34 2005-11-03 NA
## 35 2005-11-04 NA
## 36 2005-11-05 NA
## 37 2005-11-06 NA
## 38 2005-11-07 NA
## 39 2005-11-08 NA
## 40 2005-11-09 NA
## 41 2005-11-10 NA
## 42 2005-11-11 NA
## 43 2005-11-12 NA
## 44 2005-11-13 NA
## 45 2005-11-14 NA
## 46 2005-11-15 NA
## 47 2005-11-16 NA
## 48 2005-11-17 NA
## 49 2005-11-18 NA
## 50 2005-11-19 NA
## 51 2005-11-20 NA
## 52 2005-11-21 NA
## 53 2005-11-22 NA
## 54 2005-11-23 NA
## 55 2005-11-24 NA
## 56 2005-11-25 20.19373
## 57 2005-11-26 25.36564
## 58 2005-11-27 26.39569
## 59 2005-11-28 28.41766
## 60 2005-11-29 33.74083Das Ergebnis schaut gut aus, wieder sind die
Spaltenbeschriftungen nicht sehr aussagekräftig, wir
ändern diese deshalb unter Verwendung der setNames()
Funktion:
# Spaltenbeschriftung anpassen
daily_simulated_SWE = setNames(daily_simulated_SWE,c("Date","SWE"))Wir schauen kurz auf die Anzahl der Messungen und Simulationen, in
unserem Fall sollten diese beiden Zeitreihen gleich lag sein, wäre das
nicht so müssten wir einen Überlapppungszeitraum verwenden! Wir
verwenden hierzu die Funktion nrow(), welche die
Anzahl der Zeilen eines DataFrames
ausgibt:
# Anzahl der Zeilen ausgeben
nrow(daily_simulated_SWE)## [1] 273nrow(daily_observed_SWE)## [1] 273Die beiden Zeitreihen sind gleich lang (auch wenn das nicht unbedingt heisst, dass die Tage die gleichen sind), wir erstellen uns nun eine neue Liste mit den SWE Beobachtungen und Simulationen, nur der Übersichtlichkeit wegen:
# Liste mit Beobachtungen und Simulationen anlegen
scatter_input = NULL
scatter_input$Observed_SWE = daily_observed_SWE$SWE
scatter_input$Simulated_SWE = daily_simulated_SWE$SWE
# Struktur ausgeben
str(scatter_input)## List of 2
## $ Observed_SWE : num [1:273] 0 0 0 0 0 0 0 0 0 0 ...
## $ Simulated_SWE: num [1:273] NA NA NA NA NA NA NA NA NA NA ...Diese Liste übergeben wir nun an die plot() Funktion,
dabei plotten wir die Beobachtungen auf die X-Achse und die Simulationen
auf die Y-Achse. Dabei können wir Vieles anwenden, das wir eben über das
Formatieren von Plots gelernt haben:
# Scatterplot erzeugen
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, type="p", col="black", ann=FALSE)
# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")
# Grid hinzufügen
grid(lty = 'dotted', col = "gray", lwd = 1)
Unser Scatterplot schaut schon ganz gut aus, ein paar Details wollen wir aber noch ändern:
Bei einem Scatterplot muss der Wertebereich der X- und Y-Achse immer gleich sein, nur so wäre das perfekte Modell eine Linie im 45 Grad Winkel!
Wir wollen den Wertebereich auf minimal 0 umstellen, negative SWE-Werte gibt es nicht und somit brauchen wir auch keinen Bereich < 0
Die Kreissymbole wollen wir größer machen und auch füllen
Den Wertebereich können wir mit den bereits bekannten Parameter
xlim bzw. ylim anpassen, damit ist sicher
gestellt, dass X- und Y-Achse gleich skaliert sind.
Merke:
Wenn wir später unseren Plot speichern, stellen wir durch Wahl einer
Höhe und Breite sicher, dass unser Plot bei gleicher X- und Y-Skalierung
dann auch quadratisch ist, also dass das Verhältnis 1:1 sichergestellt
ist. Wenn wir in RStudio selbst plotten, ist das Grafikfenster “Plots”
je nach Rechner/Bildschirm unterschiedlich groß (wir können dieses ja
sogar größer und kleiner ziehen, was den Plot neu skaliert). Wenn wir
auch hier in der Polt-Ansicht in jedem Fall wollen, dass unser Plot in
einem Achsenverhältnis 1:1 dargestellt wird, können wir das durch das
Setzen von par(pty="s") for dem Plot Befehl tun. Alternativ
gibt es einen Parameter asp=, wenn wir diesen im Plotaufruf
auf 1 setzen, wird das Verhältnis zwischen X- und Y-Achse auf 1 gesetzt,
aber dieser Parameter ignoriert unsere xlim bzw.
ylim Einstellungen und ist somit hier nicht hilfreich.
Die Art des Symbols können wir über den Parameter
pch einstellen. Folgende Optionen sind hier möglich (wollt
Ihr diese Symbole als Grafik sehen oder mehr über den Umgang mit diesen
erfahren besucht: https://r-charts.com/base-r/pch-symbols/):
pch = 0,squarepch = 1,circlepch = 2,triangle point uppch = 3,pluspch = 4,crosspch = 5,diamondpch = 6,triangle point downpch = 7,square crosspch = 8,starpch = 9,diamond pluspch = 10,circle pluspch = 11,triangles up and downpch = 12,square pluspch = 13,circle crosspch = 14,square and triangle downpch = 15, filled squarepch = 16, filled circlepch = 17, filled triangle point-uppch = 18, filled diamondpch = 19, solid circlepch = 20,bullet (smaller circle)pch = 21, filled circle bluepch = 22, filled square bluepch = 23, filled diamond bluepch = 24, filled triangle point-up bluepch = 25, filled triangle point down blue
Die Größe des Symbols kann über cex
eingestellt werden, wobei der Wert 1 den Standard repräsentiert, größere
(z.B. 1.5) oder kleinere (z.B. 0.5) Werte vergrößern/verkleinern
entsprechend.
# Scatterplot erzeugen
par(pty="s")
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, xlim=c(0,600), ylim=c(0,600), type="p", pch = 19, cex=1.2, col="black",ann=FALSE)
# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")
# Grid hinzufügen
grid(nx = NULL, ny = NULL, lty = 'dotted', col = "gray", lwd = 1)
Um den Zusammenhang zwischen Simulationen und Beobachtungen statistisch zu beschreiben, braucht es (hier) eine lineare Regressionsfunktion. Diese wollen wir uns im Folgenden näher ansehen.
8 - Lineare Regression
Die lineare Regression versucht einen Y-Wert in linearer Abhängigkeit eines X-Wertes zu berechnen.
Heraus kommt eine Regressionsfunktion:
y = Achsenabschnitt + Steigung * x = b + m * x
Solche Regressionsfunktionen kommen zum Beispiel bei Streudiagrammen (wie hier) oder bei Trendanalysen in Klimawandelstudien zum Einsatz. In Letzteren versucht man damit die Veränderung einer abhängigen Variablen (z.B. der Temperatur, Y-Achse) in Abhängigkeit einer unabhängigen Variablen (den Jahren, X-Achse) zu untersuchen.
Um die lineare Beziehung zwischen Simulationen und
Beobachtungen herzuleiten verwenden wir die Funktion lm()
(linear model) und übergeben dieser zunächst die Y-Werte (die
Simulationen) und dann die X-Werte (die Beobachtungen). Die Ergebnisse
schreiben wir in die Variable fit_model:
# Lineare Regression berechnen
fit_model = lm(scatter_input$Simulated_SWE ~ scatter_input$Observed_SWE)Durch diesen Aufruf haben wir eine Variable fit_model
definiert, die alle Information zur berechneten Regression enthält. Die
Koeffizienten (Achsenabschnitt = Intercept) und die Steigung (hier:
scatter_input$Observed_SWE) kann man sich über Aufruf von
coefficients() und Angabe der Regressionsfunktion
fit_model ausgeben lassen:
# Koeffizieneten ausgeben
coefficients(fit_model)## (Intercept) scatter_input$Observed_SWE
## 20.759101 1.121946Was heisst das nun? Das bedeutet, ausgehend von einem Wert von 20 (mm SWE) an der Y-Achse steigt unsere Regressionsgerade mit einer Steigung von 1.12 mit jedem weiteren beobachteten mm SWE! Aber schauen wir uns das später gleich im Plot genauer an…
Will man sich gleich die gesamte Zusammenfassung der Linearen
Regression ansehen, verwendet man die Funktion
summary() unter Angabe der Regressionsfunktion
fit_model.
Ausgeben werden dann:
- unter
Calldie verwendete Funktion - die Statistik der Abweichungen (Residuen) unter
Residuals - die Regressionskoeffizienten (Achsenabschnitt und Steigung) zusammen
mit den Ergebnissen des Signifikanztests unter
Coefficients - der Standardfehler, das R2 und die F-Statistik.
Aber was heisst Signifikanz eigentlich: Der hier verwendete Signifikanztest testet die Regressionskoeffizienten gegen “0”. Er testet also, mit welcher Wahrscheinlichkeit ich richtig liege dass z.B. die Steigung ungleich “0” ist, also ein deutlicher Zusammenhang besteht.
Wir schauen uns nun die Zusammenfassung der Regression einmal an:
# Zusammenfassung ausgeben
summary(fit_model) ##
## Call:
## lm(formula = scatter_input$Simulated_SWE ~ scatter_input$Observed_SWE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -47.062 -27.555 -9.564 8.808 106.406
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.75910 7.18493 2.889 0.00441 **
## scatter_input$Observed_SWE 1.12195 0.02784 40.299 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 39.73 on 157 degrees of freedom
## (114 observations deleted due to missingness)
## Multiple R-squared: 0.9118, Adjusted R-squared: 0.9113
## F-statistic: 1624 on 1 and 157 DF, p-value: < 2.2e-16Wir sehen, dass ein signifikanter Zusammenhang besteht, das Signifikanzniveau ist mit mind. zwei Sternen angegeben, also min. 0.001. Wir gehen also mit einer Irrtumswahrscheinlichkeit von 0.1 % fälschlicherweise von einem Zusammenhang der beiden Größen aus, wir liegen also mit einer Wahrscheinlichkeit von 99.9% richtig wenn wir einen Zusammenhang annehmen!
Wir wollen nun die Regressionslinie auch im Plot
einfügen und ansehen. Dazu plotten nochmal unser Streudiagramm
und verwenden die Funktion abline() um die Regressionslinie
hinzuzufügen. Wir geben dieser einen gestrichelten Linientyp
lty='dashed', eine Linienstärke von 3 lwd=3
und die Farbe rot col='red'.
# Scatterplot erzeugen
par(pty="s")
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, xlim=c(0,600), ylim=c(0,600), type="p", pch = 19, cex=1.2, col="black", ann=FALSE)
# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")
# Grid hinzufügen
grid(lty = 'dotted', col = "gray", lwd = 1)
# Regressionslinie hinzufügen
abline(fit_model, lty = 'dashed', col = 'red', lwd = 3)
Zu einer Regressionsgerade gehört aber immer die zugehörige Regressionsfunktion, und am besten auch das Bestimmtheitsmaß (R2), beides wollen wir nun als Text hinzufügen.
Dazu müssen wir zunächst durch Zuhilfenahme der Funktion
round(Variable,Dezimalstellen) die in der Variable
coefficients gespeicherten Regressionskoeffizienten sowie
das unter r.squared zu findende Bestimmtheitsmaß
auf zwei Stellen hinter dem Komma runden.
Danach basteln wir uns eine Textvariable
regress_function in die wir mit der Funktion
paste() nacheinander Textbausteine und Variablen
einbauen.
Die Funktion text() erlaubt dann das Einfügen
der generierte Zeichenfolge in das Diagramm.
Nun das Ganze Schritt für Schritt:
Zunächst runden wir alle Variablen, stellen den Text zusammen und geben diesen Testweise ohne Diagramm aus:
# Parameter der Regressionsgleichung auf 2 Dezimalstellen runden
slope = round(summary(fit_model)$coefficients[2], 2)
intercept = round(summary(fit_model)$coefficients[1], 2)
rsquared = round(summary(fit_model)$r.squared, 2)
# Zeichenfolge mit Regressionsgleichung definieren
regress_function = paste("y = ", intercept, " + ", slope, " * x",sep='')
# Zeichenfolge mit Bestimmtheitsmaß definieren
rsquared_text = paste("R2 = ", rsquared, sep='')
# Text ausgeben
print(regress_function)## [1] "y = 20.76 + 1.12 * x"print(rsquared_text)## [1] "R2 = 0.91"Nun bauen wir diesen Text in unser Diagramm ein. Dies tun wir unter
Verwendung der text() Funktion, diese benötigt die
X-Koordinaten im Plot, die
Y-Koordinaten im Plot und den Text in
Form von:
text(X-Koordimnate,Y-Koordinate,MeinText)
Durch Hinzufügen des Parameters pos können wir noch die
Position definieren, d.h. ob der Text unter
1, links 2, über 3 oder rechts
4 von den angegeben Koordinaten platziert werden soll. Wir
vewenden den Wert 4 (rechts).
Wollt Ihr mehr über die text-Funktion erfahren, besucht: https://www.rdocumentation.org/packages/graphics/versions/3.6.2/topics/text
Die X- und Y-Koordinaten sind dabei über dei Werte der Achsen zu skalieren, sie stellen also Werte in der Plotfläche dar!
Nach Abgleich mit unseren Achsenwerten setzen wir die
Regressionsfunktion auf die X-koordinate 0 und die
Y-Koordinate 550. Den Text mit dem R2 setzen wir knapp
darunter auf die X-koordinate 0 und die Y-Koordinate
500. Bis auf die text-Funktion ganz unten verwenden wir
unseren vorherigen plot-Aufruf!
Tipp:
Während die Angabe der Koordinaten nach Sichtung des Achsenbereichs
gut über absolute Werte klappt, geht das natürlich nicht mehr wenn sich
der Wertebereich ändert! Man kann mit etwas Aufwand die Koordinaten
relativ zum Wertebereich setzen, natürlich nur wenn man diesen vorher
über xlim bzw. ylim definiert, also kennt.
Wenn ich meinen Text z.B. immer auf 95% des Wertebereiches (der Y-Achse)
setzen will kann ich das so tun:
ypos = ymin + (0.95 * (ymax - ymin))
Analog funktiomiert das auch mit der X-Achse!
# Scatterplot erzeugen
par(pty="s")
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, xlim=c(0,600), ylim=c(0,600), type="p", pch = 19, cex=1.2, col="black", ann=FALSE)
# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")
# Grid hinzufügen
grid(nx = NULL, ny = NULL, lty = 'dotted', col = "gray", lwd = 1)
# Regressionslinie hinzufügen
abline(fit_model, lty = 'dashed', col = 'red', lwd = 3)
# Regressionsfunktion als Text hinzufügen
text(0, 550, regress_function, pos=4)
# Bestimmtheitsmaß als Text hinzufügen
text(0, 520, rsquared_text, pos=4)
Das hat gut geklappt! Allerdings ist das Bestimmtheitsmaß allein nur wenig geeignet um Modellergebnisse zu bewerten, wir wollen im Folgenden deshalb weitere Effizienzkriterien kennenlernen!
9 - Effizienzkriterien
Effizienzkriterien können herangezogen werden, um Modellergebnisse zu bewerten. Folgende Kriterien sollten wir kennen:
Bestimmtheitsmaß
R2, Wertebereich: 0 <–> 1:R2 (das Quadrat des Korrelationskoeffizienten R) sagt uns welcher Teil der Gesamtvarianz durch die abhängige Variable (im Plot: Y-Achse) erklärt wird. In unserem Fall erklärt das Modell 91% der Gesamtvarianz. R2 hat aber große Nachteile, da es nur die Kovarianz betrachtet, aber nie die absoluten Differenzen zwischen Beobachtung und Modell. So kann das R2 1 sein, wenn ein Modell immer 1000 mm SWE über den Beobachtungen liegt, da beide auch dann perfekt kovarieren, mit einer guten Modellperformanz hat das natürlich dann nichts zu tun.
Nash-Suttcliffe Koeffizient
NSE, Wertebereich: minus undendlich <–> 1:Der Nash-Sutcliffe Koeffizient (auch Nash-Sutcliffe Model Efficiency, NSE) betrachtet auch die Differenzen zwischen Beobachtung und Modellergebnis und ist somit kritischer und in der Hydrologie ein Standard Kriterium. Das Perfekte Modell hätte einen NSE Wert von 1, ab einem Wert von 0 ist der Mittelwert eine bessere Vorhersage als das Modell.
PBIAS, Wertebereich: minus unendlich <–> plus unendlich:Der PBIAS gibt den Grad der Über-/Unterschätzung durch das Modell an. Ein PBIAS von -20 beschreibt eine Unterschätzung von 20%, ein PBIAS von +20 beschreibt eine Überschätzung von 20%. Das perfekte Modell liefert einen PBIAS Wert von 0.
Das waren nur einige wichtige Vertreter, für einen umfassenden Überblick zu verschiedenen Effizienzkriterien in der Hydrologie, schaut Euch folgende Artikel an:
https://adgeo.copernicus.org/articles/5/89/2005/
https://swat.tamu.edu/media/1312/moriasimodeleval.pdf
Um NSE zu berechnen könnt Ihr die Bibliothek
topmodel verwenden. Die Formel zur Berechnung von NSE
lautet:
\[
NSE=1 - [\sum_{i = 1}^{n}{(obs_i-sim_i)^2}: \sum_{i =
1}^{n}{(obs_i-obs_m)^2}]
\] Dabei stehen obs und sim für die
Beobachtungen und Simulationen, i und m für
die Werte eines Zeitschrittes bzw. das Mittel aller Werte.
Verwendet doch gleich einmal die enthaltene Funktion
NSeff(Qobs,Qsim) um NSE zu berechnen und setzt
NSE noch unter dem R2 in Euren Plot ein!
Zur Berechnung des PBIAS habe ich Euch eine kleine
Funktion gebaut, die mit PBIAS(Qobs,Qsim) gerufen werden
kann. Sie tut nichts anderes, als für jeden Zeitschritt die Abweichung
des Simulationswertes vom Beobachtungswert zu berechnen und diese
Abweichungen dann aufzusummieren. Die Summe der Abweichungen wird dann
geteilt durch die Summe aller Messwerte, mit 100 Multipliziert ergibt
sich die durchschnittliche prozentuale Abweichung:
\[ PBIAS=100 * [\sum_{i = 1}^{n}{(sim_i-obs_i)} : \sum_{i = 1}^{n}{(obs_i)}] \]
Die unten aufgeführte Funktion müsst Ihr in Eurem Skript einbauen, nachdem man die Definition einmal ausgeführt hat, kann man die Funktion immer wieder aufrufen.
# define function
PBIAS = function(Qobs,Qsim) {
# get indices with valid values
valid_obs=which(!is.na(Qobs))
valid_sim=which(!is.na(Qsim))
valid_values=intersect(valid_obs,valid_sim)
# get bias
bias=Qsim[valid_values]-Qobs[valid_values]
sum_bias=sum(bias)
sum_obs=sum(Qobs[valid_values])
pbias=100*sum_bias/sum_obs
# return result
return(pbias)
}
# use funtion to derive PBIAS
pbias=PBIAS(scatter_input$Observed_SWE,scatter_input$Simulated_SWE)
pbias## [1] 21.14462Wir überschätzen das beobachtete SWE also im Schnitt um 21%.
Verwendet nun diese Funktion um Euren PBIAS zu berechnen
und setzt den berechneten Wert noch unter dem NSE in Euren
Plot ein!
Eine grobe Richtlinie zur Bewertung von Modellergebnissen über NSE und PBIAS geben Morasi et al. (2007):
10 - Plots speichern
Für das Speichern von Plots gibt es verschiedene Möglichkeiten.
Wir können das händisch in RStudio tun, indem wir in der Plotansicht auf
Export->Save as ImagegehenWenn wir im Jupyter arbeiten, könnt Ihr die generierten Grafiken auch mit
Rechtsklick --> Speichern unter ...im Browser speichernAlternativ (und für Skripte meist besser weil automatisert möglich) können wir den Plot im Code selbst durch den
png()Befehl (vor dem Plotten einzufügen) gefolgt vondev.off()(nach dem Plotten einzufügen) speichern
Option 3) probieren wir gleich mal anhand eines Plots aus, wir
verwenden dafür nochmal den oben generierten Scatterplot, geben dabei
der png() Funktion auch noch die Breite und Höhe
unseres pngs mit, da die Qualität (Auflösung) standardmäßig
redzuiert ist. Beachtet dass wir durch eine geeignete Wahl der Breite
und Höhe dafür sorgen müssen, dass unser Achsenverhältnis nicht
verzerrrt ist, wir wählen eine leicht erhöhte Höhe um dem Titel
Platz zu geben:
# Datei aufmachen
png(file = "Scatterplot.png", width=500, height=520)
# Scatterplot erzeugen
par(pty="s")
plot(scatter_input$Observed_SWE, scatter_input$Simulated_SWE, xlim=c(0,600), ylim=c(0,600), type="p", pch = 19, cex=1.2, col="black", ann=FALSE)
# Diagrammtitel hinzufügen
title(main="Simuliertes vs. beobachtetes Schneewasseräquivalent \nStation Col de Porte 2005/2006")
# Beschriftung für X- und Y-Achse hinzufügen
title(xlab="Beobachtetes SWE (mm)")
title(ylab="Simuliertes SWE (mm)")
# Grid hinzufügen
grid(nx = NULL, ny = NULL, lty = 'dotted', col = "gray", lwd = 1)
# Regressionslinie hinzufügen
abline(fit_model, lty = 'dashed', col = 'red', lwd = 3)
# Regressionsfunktion als Text hinzufügen
text(0, 550, regress_function, pos=4)
# Bestimmtheitsmaß als Text hinzufügen
text(0, 520, rsquared_text, pos=4)
# NSE als Text hinzufügen
library(topmodel)
nse_text=paste("NSE = ", round(NSeff(scatter_input$Observed_SWE,scatter_input$Simulated_SWE),2), sep='')
text(0, 490, nse_text, pos=4)
# PBIAS als Text hinzufügen
pbias_text=paste("PBIAS = ", round(PBIAS(scatter_input$Observed_SWE,scatter_input$Simulated_SWE),2), sep='')
text(0, 460, pbias_text, pos=4)
# Device schliessen
dev.off()## quartz_off_screen
## 2Die so erzeugte Datei Scatterplot.png liegt nun im
Arbeitsverzeichnis und sollte etwa so aussehen:
